Hexagonal convection patterns in atomistically simulated fluids 
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Molecular dynamics simulation has been used to model pattern formation in three-dimensional 
Rayleigh-Benard convection at the discrete-particle level. Two examples are considered, one in 
which an almost perfect array of hexagonally-shaped convection rolls appears, the other a much 
narrower system that forms a set of linear rolls; both pattern types are familiar from experiment. 
The nature of the flow within the convection cells and quantitative aspects of the development of 
the hexagonal planform based on automated polygon subdivision are analyzed. Despite the micro- 
scopic scale of the system, relatively large simulations with several million particles and integration 
timesteps are involved. 

PACS numbers: 47.54.-r, 02.70.Ns, 47.27.T- 
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The spontaneous emergence of flow patterns in exter- 
nally driven fluids is one of the more fascinating aspects 
of fluid dynamics, and a great deal of experimental and 
theoretical effort has been spent studying such phenom- 
ena over the past century. Theoretical fluid dynamics has 
successfully ignored the underlying atomicity, describing 
the dynamics in terms of nonlinear partial differential 
equations for continuous fields; analysis of flow instability 
within this framework demands approximations whose 
consequences are not always readily assessed. Ideally, 
one would like to understand flow instability in terms of 
the behavior at the microscopic level by looking beyond 
the continuum representation to address the dynamics of 
the molecules themselves, albeit at considerable compu- 
tational cost; given sufficiently many molecules followed 
for an adequate time interval, one can reasonably expect 
to see the same kind of behavior (initially in a qualita- 
tively correct form, then, with even larger systems, also 
quantitatively) . The ability to simultaneously encompass 
the dynamics both at the molecular level and the scale 
at which the cooperativity underlying structured flow be- 
comes manifest is an important step in trying to break 
free of the limitations of traditional continuum fluid dy- 
namics. 

Accomplishing this goal calls for molecular dynamics 
(MD) simulation, a technique in broad general use, but 
less so for fluid dynamics (quite likely due to Avogadro's 
tyranny - an apparent need for extremely large systems 
to accommodate the phenomena). A very early MD effort 
related to bridging the gap between atomistic dynam- 
ics and fluid behavior at the continuum level was under- 
standing long-time effects due to correlations in atomic 
trajectories [J, but it was only much later that actual 
MD simulation of complex flow was attempted. Initial 
efforts were confined to two dimensions, starting with the 
wakes of obstructed flows 0,01 and followed by Rayleigh- 
Benard (RB) convection Jill H Sll @; the study of 
Taylor-Couette flow 0, demonstrated that complex 
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three-dimensional flow could also be dealt with and, sub- 
sequently, the Rayleigh-Taylor instability was also mod- 
eled 0]. Other instances of MD contributing to fluid 
dynamics include the moving interface between immis- 
cible liquids undergoing Poiseuille flow 0] and droplet 
breakup in liquid jets |l4| . The present study describes 
the application of the MD approach to the full three- 
dimensional RB problem, where the behavior is poten- 
tially far richer than in two dimensions, and where simu- 
lation is more closely related to experiment; the ability to 
tackle increasingly demanding MD simulations, both in 
terms of size and duration, reflects the availability of in- 
creasingly powerful and affordable computing resources. 

The RB phenomenon - the rich variety of flow pat- 
terns produced by convection in a fluid layer heated 
from below - continues to be a problem of great in- 
terest 0, 0, 0, 0, 0] . The principal dimensionless 
quantity governing the behavior is the Rayleigh number, 
Ra = agL^AT/vK, where a, v and k are the thermal ex- 
pansion coefficient, kinematic viscosity and thermal dif- 
fusivity, L z is the layer height, g the gravitational accel- 
eration, and AT = T\ — To is the temperature difference 
between the lower hot wall T± and the upper cold wall 
To; a critical value Ra c marks the onset of convection, 
where buoyancy overcomes viscous drag, and convection 
replaces conduction as the preferred thermal transport 
mechanism. 

Theory is simplified by the Boussinesq approximation, 
in which it is assumed the only temperature-dependent 
fluid property is the density. Ra c can be computed 
for various kinds of boundary conditions |15| . as can 
the wavelength A of the convection pattern at critical- 
ity, but not the actual planform (e.g., rolls, spirals and 
hexagons). Roll width (typically ~ L z ) represents a com- 
promise; wide rolls reduce both viscous drag and diffusive 
heat transfer between ascending hot and descending cold 
streams, whereas narrow rolls reduce shear losses close 
to nonslip walls. Hexagonal patterns were once associ- 
ated with non-Boussinesq fluids 0, |2(| (and are better 
known in convection driven by surface tension) ; the flow 
direction at the cell center is determined by the temper- 
ature dependence of v and differs for gases and liquids. 



2 




FIG. 1: (Color online) Final streamline plot for the large rect- 
angular system #A showing the hexagonal cell pattern; the 
streamlines are color coded to indicate temperature variation 
(ranging from red for hot to blue for cold). 



Hexagons are now known to occur in (essentially) Boussi- 
nesq systems above Ra c 21, 22], even forming coexisting 
regions with opposite flows at the cell centers, a phe- 
nomenon not attributable to variable v. 

The MD simulation [2^ considers a fluid of soft-sphere 
atoms with a short-ranged, repulsive interaction u(r) = 
4e[(cr/r) 12 - (a/r) 6 ], with a cutoff at r c = 2 1 / 6 cr. In the 
dimensionless MD units used subsequently, a and e de- 
termine length and energy, while temperature is defined 
by setting ks = 1; in the case of Argon, a = 3.4 A, 
e/ks = 120K, and the unit of time is 2.16ps. The top 
and bottom thermal walls of the system are each formed 
from a layer of fixed atoms arranged as a lattice, spaced 
to ensure roughness and impenetrability; lateral bound- 
aries are periodic. The temperature gradient is produced 
by rescaling the velocities of those atoms adjacent (within 
a range of 1.3) to the thermal walls; rescaling occurs ev- 
ery 20 timesteps to allow the effect to propagate without 
unduly affecting the dynamics, and is applied to atoms 
moving away from the walls. In the initial state, atoms 
are placed on the sites of a regular grid, with random ve- 
locities corresponding to a uniform vertical temperature 
gradient. The integration timestep is 0.004. In view of 
the extensive computations required, parallel processing 
based on spatial decomposition is used. 

Flow analysis entails spatial coarse-graining, in which 
a grid subdivision is applied to the region and the mean 
properties (velocity, temperature, density) for each grid 
cell are evaluated for its occupant atoms. This is re- 
peated every 20 timesteps and averages over 100 succes- 
sive measurements recorded, with typical grid sizes of 5-6 
horizontally and 3 vertically (rj 30 atoms per grid cell); 
this represents a compromise between capturing the fine 
spatial and temporal detail of the developing convection 
patterns and reducing noise due to thermal fluctuations. 

Two of a series of runs are described in detail below; 
the eventual outcomes of both are organized convection 
cell patterns, a hexagonal array in the case of system 




FIG. 2: (Color online) Early temperature distribution near 
the lower wall. 



#A and a set of linear rolls for #B. System #A contains 
N a = 3 507 170 atoms (of which N w = 447458 are fixed 
in the walls) and is run for a total of 3.19 x 10 6 timesteps 
(= 2.8 x 10~ 8 s). The region size is L x x L y x L z = 
521 x 451 x 35.3; despite the smallness of L z (a mere 
120 A), the layer thickness exceeds the value 25 used in 
the Taylor-Couette study 0, 0] . Since the rate of pat- 
tern change slows as features become larger, a reflection 
of the underlying horizontal diffusion processes, the long 
run is needed to ensure a final steady state is reached. 
The smaller system #B has L y = 65 (L x and L z are un- 
changed) , resulting in a more modest N a = 506 696 (with 
N w = 64328), and a run length of 0.48 x 10 6 timesteps. In 
both cases, T = 1, Ti = 10; the value g = 3AT/2L Z en- 
sures equality of the potential and kinetic energy changes 
across the layer and avoids excessive density variation. 
Large gradients and fields are the norm for MD flow stud- 
ies that must deal with relatively small regions, here com- 
pensating for the small L z to ensure an adequate value 
of Ra; AT/Tq is large compared with experiment (where 
AT/To w 0.01 or less), and g = 2.7 x 10 12 g earth (the 
variation in gravitational potential across the layer, gL z , 
is only 8 times that of The total computation time 

for system #A amounts to approximately 6 cpu-months 
on a 2 GHz Intel Pentium processor, with the actual run 
using an 8-cpu cluster. 

The final state of system #A is shown in Fig. ^ as 
a streamline plot that reveals the full three-dimensional 
flow structure; the curved tracks (color coded to show 
temperature variation) are derived from the coarse- 
grained grid (size 96 x 82 x 12). The most prominent 
feature is the array of hexagonal convection cells extend- 
ing across the full depth of the layer, with cold fluid de- 
scending at the cell centers. The cell array is free to 
orient itself to accommodate a preferred A for the partic- 
ular size and shape of the periodic container; animated 
three-dimensional visualization allows following the pat- 
tern development as cells grow and merge (not shown). 

Figure [3 shows the horizontal temperature distribu- 
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FIG. 3: (Color online) Final temperature distribution with 
Voronoi polygon subdivision superimposed. 
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FIG. 4: Graphs showing time-dependence (dimensionless MD 
units) of the cell count and the mean radial range (typical 
fluctuations are included for the latter). 

tion at a distance « L z /5 from the hot wall early in the 
run, after the onset of convection but before a regular 
cell pattern develops. The distribution at the end of the 
run appears in Fig. |3 with the periodic cell array clearly 
visible. Superimposed on this image is a polygon subdi- 
vision produced by an automated Voronoi analysis of a 
slice through the flow structure (limited system size can- 
not accommodate the large aspect ratios L x /L z > 100 
used experimentally, as opposed to 15 here, that yield 
extended patterns amenable to Fourier and other forms 
of analysis [24|); cell centers marked by white spots are 
located at each of the local temperature minima (or at 
the downward flow maxima) and the plane subdivided 
into polygons in which the region nearer to a particu- 
lar center than to any other is assigned to its polygon. 
The resulting polygon set provides a good fit to the cell 
boundaries where hot upward flow is strongest (the fit 
quality improves as the run progresses), and is used in 
the quantitative analysis below. 
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FIG. 5: Flow field azimuthally averaged over all cells; the 
axes measure the radial distance from the cell center and the 
height above the hot wall (dimensionless MD units). 



The time-dependence of two prominent features of the 
polygon arrays is shown in Fig. The polygon count 
falls from an initial value of w 30 to the eventual 8 at 
time t ss 5000, it then increases to 9 and finally drops 
back to 8 at t ss 8000 as cells split and merge; this be- 
havior is an indication of the slow relaxation processes 
that the simulations must be able to cover. The mean 
deviation of the polygon shape from regularity, expressed 
as the range of radial distances of the nearest and fur- 
thest vertices of each polygon, shows a more gradual con- 
vergence to the final value of 4.5, which is just 5% of 
the mean radius. Other planform properties include the 
mean distance from the center to the nearest point on 
an edge, \hex/2, whose final value is 93 ± 1, the mean 
polygon area that is inversely proportional to the count, 
and the area variation that falls to a minimum of ± 2% 
of the mean in the final state. The average number of 
vertices per polygon must be 6; the spread of values falls 
to zero as the final state is neared. These results quantify 
what is observed directly, namely an initially disordered 
set of small convection cells that grow, merge and rear- 
range, culminating in what appears to be a stable array 
of hexagonal cells; the final wavelength Xtex = 5.2 7 Lx , a 
value not much larger than the experimental 4 L z |2l| (a 
lower Xhex is reported in |22^): the characteristic size of 
the cells present in the early, irregular pattern is about 
half this value (Fig.|2J). 

Details of the flow profile inside a single convection cell 
appear in Fig. [5] The arrow plot shows the radial and 
vertical components of the flow field within the cells, av- 
eraged over all cells of the final state, and is evaluated 
by projecting the flow velocity onto a series of vertical 
planes arranged as radial spokes around each cell center, 
with only grid cells that intersect the planes contribut- 
ing. The maximum arrow length corresponds to a speed 
of 0.52. While average descending flow appears faster 
than ascending, this is consistent with the toroidal roll 
shape that allows more space for the latter. The nonslip 
nature of the flow at the lower wall is partly masked by 
coarse-graining, but flow rate is seen to decrease near the 
wall (the effectiveness of the nonslip mechanism is bet- 
ter appreciated in small systems run interactively). Lo- 
cal temperature (not shown) increases radially outwards 
from the cell core and drops with height; its estimation 
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FIG. 6: (Color online) Final streamline plot for system #B 
showing the linear roll pattern. 

is problematic because the contribution of nonuniform 
bulk flow cannot be fully removed, and thermalization in- 
cludes atoms interacting with the walls, so the measured 
range is 0.84-11.4 (rather than 1-10), but the gradual 
variation across the roll profile and the absence of high 
gradients near the walls (except where the cool downflow 
initially encounters the hot wall) are evidence of efficient 
heat exchange between fluid and walls. Local density 
varies from 0.58 at the cell core to 0.38 at the periphery 
(the overall mean is 0.4). 

The final state of run =#=B consists of a set of 6 pairs of 
counter-rotating rolls. Fig. shows the streamline plot; 
the wavelength X ro ii (2x roll diameter) = 87 = 2.46 L z 
(note that 7 roll pairs would have produced a value close 
to the theoretical X ro u = 2.016 L z 0]). The wavelength 
ratio from runs $A and #B is Xhex/X ro ii — 2.14; the 
experimental ratio for coexisting hexagons and rolls is 
1.2—1.3 |2l|. At an early stage in run #B there is an 
unsuccessful attempt to form a pair of longitudinal rolls. 

When T\ is lowered to 6 the roll pattern in system #B 
is unchanged, but at T\ = 4 the results are noisy and roll 
structure is not readily discerned, though it is still present 
in weakened form. Strong finite-size effects due to small 
L z erase any hint of a sharp transition between conduct- 
ing and convecting states. The fluid has its own intrinsic 
thermodynamic and transport properties determined by 



the soft-sphere model and, since AT /To is large, does not 
satisfy the Boussinesq condition; estimating Ra based on 
Enskog theory and a simple equation of state 0, 13 is 
inappropriate here due to the variation of v and k across 
the layer. 

These simulations (and work in progress) provide 
strong confirmation that MD is capable of reproducing 
known fluid dynamical behavior; the unique aspect of 
this approach is that, while following the dynamics at the 
atomic scale where the underlying causes of flow insta- 
bility and self-organization are to be found, it simultane- 
ously embraces the continuum scale where fluid behavior 
is manifest and scant evidence of atomism remains. In- 
creasing computer performance (parallel computers are 
ideal for simulations of this kind) allows ever-larger sys- 
tems to be modeled, although this growth is tempered 
by the fact that computation time increases not only 
in proportion to system size but also due to slow diffu- 
sive processes that govern structure development. Larger 
systems reduce the necessarily high thermal and veloc- 
ity gradients and help eliminate spurious finite-size ef- 
fects, allowing closer correspondence with experiment. 
The outcome of the present MD study has obvious im- 
plications for future MD exploration of fluid dynamical 
phenomena. 
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